### To work out population viability, and mean dispersal evolving across a range of habitat, lambda, and d.cost values
### All values assessed over 500 generations on a 10x10 lattice with K=40

source(file="/home1/30/jc227089/evo-dispersal/SRE/Stable range functions.R", echo=F)
setwd("/home1/30/jc227089/SRE/equil/ran_outs")
# to collect data
collect<-function(popmatrix, ...){
	if (length(popmatrix[,1])==0){
		popsize<-0
		mean.disp<-NA
		sd.disp<-NA	
	}	
	else {
		popsize<-length(popmatrix[,1])
		mean.disp<-mean(popmatrix[,"P"])
		sd.disp<-sd(popmatrix[,"P"])	
	}
	out<-cbind(..., popsize, mean.disp, sd.disp)
	out
}

### initialise
args=(commandArgs(TRUE))

#evaluate the arguments
# input arguments will be number of replicates runs (rr), file.ID, habitat value (hv) and dispersal cost (dc)
for(i in 1:length(args)) {
	 eval(parse(text=args[[i]]))
}

nbrs<-neighbours.init(10, 10) #neighbour matrix
data<-c()
for (r in 1:10){ # lambda 1:10
	for (rep in 1:rr){ #reps
		pop<-init.inds(n=10*10*20, 10, 10) # half K individuals on 10x10 lattice
		hab<-init.hab(10,10,0,hv) # homogenous habitat with surv value = hv
		for (g in 1:500){ #500 generations
			if (length(pop[,1])==0) break
			pop<-repro.disp(pop, hab, 10, 10, K=40, lambda=r, d.cost=dc, n.list=nbrs, mutn=0.05)
		}
		out<-collect(pop, r, rep, g, hv, dc) #collect data
		data<-rbind(data, out)
	}	
}

save(data, file=paste("SRE-equil", file.ID, ".RData", sep=""))